Minimization of a univariate polynomial

Contributed by: Benoît Legat

using DynamicPolynomials
using SumOfSquares
import CSDP

Consider the problem of finding both the minimum value of p = x^4 - 4x^3 - 2x^2 + 12x + 3 as well as its minimizers.

We can use SumOfSquares.jl to find such these values as follows. We first define the polynomial using DynamicPolynomials.

@polyvar x
p = x^4 - 4x^3 - 2x^2 + 12x + 3

\[ x^{4} - 4x^{3} - 2x^{2} + 12x + 3 \]

Secondly, we create a Sum-of-Squares program searching for the maximal lower bound σ of the polynomial.

model = SOSModel(CSDP.Optimizer)
@variable(model, σ)
@constraint(model, cref, p >= σ)
@objective(model, Max, σ)

\[ σ \]

Thirdly, solve the program and find σ = -6 as lower bound:

optimize!(model)
solution_summary(model)

* Solver : CSDP

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "Problem solved to optimality."

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : -6.00000e+00
  Dual objective value : -6.00000e+00

* Work counters
  Solve time (sec)   : 1.49608e-03

We can look at the certificate that σ = -6 is a lower bound:

sos_dec = sos_decomposition(cref, 1e-4)

(-0.9999999995313261*x^2 + 1.9999999988281423*x + 3.000000004901299)^2

Indeed, p + 6 = (x^2 - 2x - 3)^2 so p ≥ -6.

Extraction of minimizers

We can now find the minimizers from the moment matrix:

ν = moment_matrix(cref)
ν.Q

3×3 SymMatrix{Float64}:
 22.2472   6.43654    3.12472
  6.43654  3.12472    0.0623621
  3.12472  0.0623621  1.0

This matrix is the convex combination of the moment matrices corresponding to two atomic measures at -1 and 3 which allows us to conclude that -1 and 3 are global minimizers.

η = extractatoms(ν, 1e-4)
minimizers = [η.atoms[1].center; η.atoms[2].center]

2-element Vector{Float64}:
  2.999999948889303
 -1.0000001972175558

Below are more details on what we mean by convex combination. The moment matrix of the atomic measure at the first minimizer is:

η1 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[1])), ν.basis.monomials)
η1.Q

3×3 SymMatrix{Float64}:
 81.0  27.0  9.0
 27.0   9.0  3.0
  9.0   3.0  1.0

The moment matrix of the atomic measure at the second minimizer is:

η2 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[2])), ν.basis.monomials)
η2.Q

3×3 SymMatrix{Float64}:
  1.0  -1.0   1.0
 -1.0   1.0  -1.0
  1.0  -1.0   1.0

And the moment matrix is the convex combination of both:

Q12 = η1.Q * η.atoms[1].weight + η2.Q * η.atoms[2].weight

3×3 Matrix{Float64}:
 22.2472   6.43654    3.12472
  6.43654  3.12472    0.0623623
  3.12472  0.0623623  1.0

Another way to see this (by linearity of the expectation) is that ν is the moment matrix of the convex combination of the two atomic measures.

Changing the polynomial basis

The monomial basis used by default can leave a problem quite ill-conditioned for the solver. Let's try to use another basis instead:

model = SOSModel(CSDP.Optimizer)
@variable(model, σ)
@constraint(model, cheby_cref, p >= σ, basis = ChebyshevBasisFirstKind)
@objective(model, Max, σ)
optimize!(model)
solution_summary(model)

* Solver : CSDP

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "Problem solved to optimality."

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : -6.00000e+00
  Dual objective value : -6.00000e+00

* Work counters
  Solve time (sec)   : 2.26092e-03

Although the gram matrix in the monomial basis:

g = gram_matrix(cref)
@show g.basis
g.Q

3×3 SymMatrix{Float64}:
  1.0  -2.0  -3.0
 -2.0   4.0   6.0
 -3.0   6.0   9.0

looks different from the gram matrix in the Chebyshev basis:

cheby_g = gram_matrix(cheby_cref)
@show cheby_g.basis
cheby_g.Q

3×3 SymMatrix{Float64}:
  0.25  -1.0  -1.25
 -1.0    4.0   5.0
 -1.25   5.0   6.25

they both yields the same Sum-of-Squares decomposition:

cheby_sos_dec = sos_decomposition(cheby_cref, 1e-4)

(-0.999999999927962*x^2 + 1.999999999846479*x + 3.000000000736021)^2

The gram matrix in the Chebyshev basis can be understood as follows. To express the polynomial $-x^2 + 2x + 3$ in the Chebyshev basis, we start by substituting $x$ into $\cos(\theta)$ to obtain $-\cos(\theta)^2 + 2\cos(\theta) + 3$. We now express it as a combination of $\cos(n\theta)$ for $n = 0, 1, 2$: $-(2\cos(\theta) - 1) /2 + 2 \cos(\theta) + 5/2.$ Therefore, the coefficients in the Chebyshev basis is:

cheby_coefs = [-1/2, 2, 5/2]

3-element Vector{Float64}:
 -0.5
  2.0
  2.5

We can indeed observe that we obtain the same matrix as cheby_g.Q

cheby_coefs * cheby_coefs'

3×3 Matrix{Float64}:
  0.25  -1.0  -1.25
 -1.0    4.0   5.0
 -1.25   5.0   6.25

This page was generated using Literate.jl.